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Abstract 



The purpose of this work is to test the application of the finite element 
method to quantum mechanical problems, in particular for solving the Schrodinger 
equation. 

We begin with an overview of quantum mechanics, the Schrodinger equa- 
tion and numerical techniques used to solve quantum mechanical problems. 
We note that one of the most important aspects of using the Crank-Nicolson 
method in solving the Schrodinger equation is that the numerical time step- 
ping equations are unitary, thus they inherently conserve probability, which 
is an important factor of quantum physics. 

We give an introduction to finite element analysis using the diffusion 
equation as an example. We consider three numerical time evolution meth- 
ods: the (tried and tested) Crank-Nicolson method, the continuous space- 
time method, and the discontinuous space-time method. Once a numerical 
background is established we apply these techniques to quantum mechanical 
problems: a wave packet trapped in an infinite quantum well, and a wave 
packet trapped in an infinite well with a finite barrier. 

The first point of interest is that the finite element equations associated 
with the continuous and discontinuous space-time methods are not unitary. 
We show that the explicit part of the continuous method is unstable, and the 
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implicit part is heavily damped, so both are as bad as using the forward or 
backward Euler methods. However, it is also shown that when the implicit 
and explicit parts are combined by taking their average we obtain the Crank- 
Nicolson method, which is stable and unitary. It is also shown that the 
discontinuous space-time method suffers from a small amount of damping, 
which can be controlled by the timestep size, however this comes at the 
cost of greater computation time. From this we conclude that the standard 
Galerkin space-time methods are not as good as the Crank-Nicolson method. 

On the other hand the Crank-Nicolson method also has its limitations. 
When a particle interacts with a barrier i.e. when a change of potential 
occurs, a fluctuation in probability conservation also occurs. This is shown 
to happen due to the wave packet splitting into a reflected and transmitted 
part, which increases the complexity of the wave function. It is shown that 
these fluctuations can be controlled by resolving the wave function more 
accurately, which is achieved by increasing the number of spatial elements 
used. This is further explored when we show that slight "damping" takes 
place when we model a particle in a sinusoidal potential, which is a result 
of the packet undergoing a chain of reflections and transmissions. In this 
case a much finer resolution (> 3000 spatial elements) is required in order to 
accurately model the wave function at later times. 
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Chapter 1 
Introduction 



The development of quantum mechanics in the early part of the 20th century 
led to a greater understanding of the atomic and subatomic world. On the 
nano-scale quantum theory is a more fundamental theory than Newtonian 
mechanics and classical electromagnetism as it provides a more accurate de- 
scription of many phenomena which are never observed in the macroscopic 
world. Everything from understanding the discrete nature of observable prop- 
erties such as energy and momentum to the conceptual leap of modelling 
particles as waves has all had a profound effect on the way we see and use 
the microscopic world. 

One of the biggest applications of quantum mechanics was in the mid 
20th century, as it formed the framework for understanding and develop- 
ment of semiconductor materials and devices, such as transistors, diodes, 
solar-cells, lasers, and microprocessors. However, the ongoing research in 
quantum mechanics and its application to solid-state physics has led to far 
more complex technologies, such as quantum dots and quantum wires. These 
new technologies offer important opportunities as the building blocks for the 
next generation of electronic and opto-electronic devices ranging from ultra- 
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fast optical switching to ultra-dense memories. Some recent devices which 
inherently employ such quantum technologies are quantum- well lasers, which 
have led to the development of the compact blue laser used in high storage 
optical disc players like the Playstation 3, and high performance solar-cells. 
With the huge potential of quantum wires and dots they are under immense 
active experimental and theoretical investigation. 

In order to test new theories, devices and applications it has become more 
popular from the business point of view to construct computational models 
before any physical experimentation takes place. Sometimes it is also useful 
to do practical and numerical tests side by side in order to compare results. 
So as the development of more complex quantum devices continues we will 
require more efficient and accurate numerical techniques. 

In the past models such as drift-diffusion formed the basis for simulat- 
ing semiconductor devices, but such techniques are not adequate to model 
the new breed of quantum devices where the quantum effects of a single 
electron can play a significant part in a device's operation. However, a con- 
cise quantum mechanical simulation of an entire semiconductor device is 
not feasible from the numerical point of view. So it has been stated that 
to model more complicated quantum devices we could use the fact that in 
many semiconductor devices quantum effects take place in a localized region 
(micro-structure), for example within the active zone of a quantum well laser, 
whereas the rest of the device (macro-structure) can be described by classical 
models [1]. Therefore, it would be possible to follow a strategy where we can 
couple quantum mechanical and macroscopic models (similar to multi-scale 
modelling in other areas of engineering). 



CHAPTER 1. INTRODUCTION 3 

An important aspect of quantum phenomena is the conservation of prob- 
ability. For this reason the promising time evolution numerical technique 
seems to be the Crank-Nicolson method, which is not only unconditionally 
stable but the time-stepping equations associated with it are unitary and 
thus inherently conserve probability 1 . Then for the numerical solution of the 
Schrodinger equation one can simply use a finite-difference method for space 
discretization and then apply the Crank-Nicolson method for time evolu- 
tion. In previous work, [11], more sophisticated techniques such as the finite 
element method have been used for spatial discretization. By using finite 
elements combined with high-order spatial discretisation we can then model 
the wave function of a particle more accurately. Combining this system with 
the Crank-Nicolson time evolution method we can then model the time evo- 
lution of a wave function efficiently and accurately, rather than simply using 
a simple finite difference method [12]. One step beyond this is the use of 
finite elements in space and time - known as the space-time finite element 
method. 

The most important aspect of our work will be the comparison of the 
Crank-Nicolson and finite element method to the space-time finite element 
method. We will work with simple micro-structure models such as particles 
(represented by wave-packets) in quantum wires (a one dimensional quantum 
well). The Crank-Nicolson and space-time models will then be compared for 
their efficiency, conservation of probability and accuracy. 

In Chapter 2 we will begin with an overview of quantum mechanics, the 

1 The unitary property and probability conservation will be discussed in more detail in 
the next chapter. 
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derivation of the Schrodinger equation, and some simple numerical solution 
methods. Chapter 3 will form the basis of our work on the Schrodinger equa- 
tion. We will begin with an introduction to the finite element method, then 
we will go on to deriving the element equations for the diffusion equation 2 
by first using the Crank-Nicolson method, then the continuous space-time 
method and finally the discontinuous space-time method. In Chapter 4 we 
go on to model the Schrodinger equation using the previously discussed meth- 
ods. We model simple systems such as infinite quantum wells, and a quantum 
well with a barrier. We show that the Crank-Nicolson method is by far more 
efficient and accurate to the space-time method, however we also show the 
limitations of the Crank-Nicolson method. In Chapter 5 we go on to anal- 
yse the limitations of the Crank-Nicolson method by applying it to a packet 
trapped in an infinite well with a sinusoidal potential. Finally in Chapter 6 
we give a summary of the work. 

Note on Notation 

In this work we will be dealing with a large number of matrices and opera- 
tors. To keep things consistent we will use the following notation: 



2 Due to its similar structure to the Schrodinger equation we are able to use the results 
in later chapters. 



X or x' 
X / 



N 
B 
I 



Finite element basis operator. 
Finite element derivative operator. 
Diagonal identity matrix. 
Intermediate step in matrix X. 
Final real matrix form of X. 



Chapter 2 

Quantum Mechanics 



In general, quantum physics is concerned with processes which involve dis- 
crete energies and quanta (i.e. single particles such as the photon). The mo- 
tion and behaviour of quantum processes can be described by the Schrodinger 
equation. The use of the Schrodinger equation to study quantum phenomena 
is known as Quantum Mechanics, akin to classical mechanics being the tool 
to study classical physics. In this chapter we will give a brief overview of 
quantum mechanics. Beginning with the postulates of quantum mechanics, 
we will go on to discuss the derivation of the Schrodinger equation and give 
some simple applications. We will end this chapter with a description of the 
current numerical techniques used to solve the Schrodinger equation. 

2.1 Postulates of Quantum Mechanics. 

• Postulate 1 The state of a quantum mechanical system is completely 
specified by a function ^>(x, t), which depends on the space and time 
coordinates of the particle. This function, called the wave function or 
state function, has the important property that its norm ^(x, t)i/j(x., t)dv 
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is the probability that the particle lies in the volume element dv located 
at x at time t. The wavefunction must satisfy certain mathematical 
conditions because of this probabilistic interpretation. For the case of 
a single particle, the probability of finding it somewhere is 1. So we 
have the normalization condition: 

/oo 
tp*(x,t)ip(-x,t)dv = 1 (2.1) 
-OO 

The wavefunction must also be single- valued, continuous, and finite. 

• Postulate 2 To every observable, A, in classical mechanics (e.g. energy 
and momentum) there corresponds a linear Hermitian operator, A, in 
quantum mechanics. 

• Postulate 3 In any measurement of the observable associated with op- 
erator A, the only values that will ever be observed are the eigenvalues 
a, which satisfy the eigenvalue equation: 

Ma = a*jj a (2.2) 

where ip a is the eigenfunction associated with the eigenvalue a of the 
operator A. This postulate captures the central point of quantum me- 
chanics that values of dynamical variables can be quantized. If the sys- 
tem is in an eigenstate of A with eigenvalue a, then any measurement 
of the quantity A will yield a. Although measurements must always 
yield an eigenvalue, the state does not have to be an eigenstate of A. 
An arbitrary state can be expanded in the complete set of eigenvectors 
of A {Aipi diipi) as: 

n 

^ = J2a^i (2-3) 
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In this case we only know that the measurement of A will yield one of 
the values ai with a probability |q| 2 

• Postulate 4 If a system is in a state described by a normalized wave 
function ^, then the average value of the observable corresponding to 
A is given by: 

/oo 
tp*(-x,t)Aip{x,t)dv (2.4) 
-OO 

• Postulate 5 The wavefunction of a system evolves in time according 
to the time-dependent Schrodinger equation: 

^-V(x, t) = -—VV(x, t) + ^(x)^(x, t) (2.5) 

where m is the mass of the particle, h — is the Planck constant 
and V(x) is a real function representing the potential energy of the 
system. Although the time-independent Schrodinger equation can be 
derived through elementary methods (discussed in the next section), 
the time-dependent version can not be derived so must be accepted as 
a fundamental postulate of quantum mechanics. 

2.2 The Schrodinger Equation 

In 1925 Erwin Schrodinger developed a method of quantum mechanics in- 
volving partial differential equations. This method differed to the one de- 
veloped earlier by Werner Heisenberg which employed matrices. These dif- 
ferential and matrix based methods were later shown to be mathematically 
equivalent [3]. 
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2.2.1 Time-Independent Schrodinger Equation 

One of the fundamental concepts of quantum physics is that of wave-particle 
duality: that is waves can behave like particles and particles like waves. For 
example, Einstein showed that a photon, which is considered to be a wave 
packet, has momentum just like a particle moving with the same energy, 
Appendix A. The dynamical behaviour of these quantum waves/particles can 
be described in a non-relativistic 1 manner through the use of wave mechanics. 
The single-particle three-dimensional time-dependent Schrodinger equation 
is given in Eqn. (2.5). Before we consider the full time-dependent equation, 
which must be accepted as a postulate of Quantum Mechanics, we will give 
a brief derivation of the time-independent version, which has a conceptual 
derivation linked to the wave equation. 

Derivation of the Time-Independent Schrodinger Equation 

Starting with the one-dimensional classical wave equation, 



d 2 u 



1 d 2 u 



(2.6) 



dx 2 



v 



and using separation of variables, 



u(x,t)=i/>(x)f(t) 



(2.7) 



we obtain 



(2.8) 



^Tor the relativistic description of particles and waves we require the Dirac equation for 
spin \ particles, the Klein Gordon equation for spin particles. This is all encompassed 
more generally in the study of Quantum Field Theory. 
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Then, using a standard solution of the wave equation, f(t) = e ZUJt , we obtain 

^(x) = ~1>(x) (2-9) 

This gives an ordinary differential equation describing the spatial amplitude 
of the matter wave as a function of position. This can be put in the standard 
form for the Schrodinger equation by using the fact that the energy of a 
particle is the sum of kinetic and potential parts, 

E=£ + V(x), (2.10) 

Finally, using uj = 27rz/, v = v\ and h = pX we have 

uj 2 _ 4tt V _ 4tt 2 _ 2m[E - V(x)} 

^ ~ v 2 h 2 ^' ' 

which when combined with Eqn. (2.9) gives 

d 2 2m 

—^(x) + W [E- V{x)]t/>{x) = (2.12) 

This single-particle one-dimensional equation can be extended to the case of 
three dimensions, where after rearranging it becomes 

h 2 

-— VV(x) + V(x)^(x) = £^(x) (2.13) 

Zi ill 

The solutions to this equation then represent the state function of a particle 
of mass m in a potential V(x). 

2.2.2 Time-Dependent Schrodinger Equation 

As stated in the previous section, although the time-independent Schrodinger 
equation can be derived analytically, the time-dependent Schrodinger equa- 
tion cannot be derived using such methods and is therefore generally con- 
sidered as a postulate of quantum mechanics [2]. However, we are able to 
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show that the time-dependent equation is a reasonable model of the dynamic 
evolution of a particle's states function even though it is not derivable. As 
before, using separation of variables, 

^(x,*)=^(x)/(*), 

and substituting this into Eqn. (2.5) we have 

ih df 1 f h 2 « T r/ J , , x 

7wi = ^r^ v+y(x r (x) (2 - 14) 

Now, as the left-hand side is a function of t only and the right hand side is 
a function of x only, the two sides must be equal to a constant. Assigning 
this constant as as the right-hand side clearly has dimensions of energy, 
we can then extract two ordinary differential equations: 

fit) dt ~ h [lAi>) 

and where the other is the time-independent Schrodinger equation, Eqn. (2.13). 
Simply solving Eqn. (2.15) we have 

f(t) = e~ iEt ' h (2.16) 

The energy operator, given by Eqn. (2.13), known as the Hamiltonian is 
a Hermitian operator, therefore its eigenvalues are real, so E is real. This 
means that the solutions of Eqn. (2.15) are purely oscillatory. Therefore, if 

^(x,t)=^(x)e-^/ ft , (2.17) 

then the total wave function ^(x, t) differs from ^( x ) on ly by a phase factor 
of constant magnitude. This then implies that the probability, or the norm, 
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of the particle state is time independent, 



|^(x,t)| 2 = ^(M)iKM) = e^/V(x)e-^V(x) = ^(x)^(x) = |^(x)| 2 

(2.18) 



It also implies that the expectation value for any time-independent operator 
is also time-independent, 



For this reason the states described by the wavefunction in Eqn. (2.17) are 
called stationary states. However, even though the probablity distribution 
described by ^( x ^) is stationary, the particle it describes is not. This could 
be conceptually understood by having a particle in a box. In such a case 
the particle will be moving around in the box: bouncing off the walls etc. 
However, the probability distribution of the particle within the box will be 
constant in time. Thus, if the probability is 0.5 in the middle of the box, 
and at the box edges, this implies that if we check for the particle in 100 
identical boxes we will find it in the middle in 50 of them. 

2.3 Analytical Solutions 
2.3.1 Particle in a Box 

As a simple example we consider a particle constrained to move in a single 
dimension under the influence of a potential V{x) which is zero for < x < 
L and infinite elsewhere, Fig. 2.1. Since the wavefunction is not allowed 
to become infinite, it must have a value of zero where V(x) is infinite i.e. 
^(0) = ip(L) = 0, so ip(x) is nonzero only within [0,L]. The Schrodinger 




(2.19) 
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V = oo 



V = 



V = oo 







Figure 2.1: Infinite Square well. 



equation for this simple case is 
h 2 (ftpix) 



2m dx 2 



= Eip{x) 



< x < L 



(2.20) 



Solving this and applying the normalization condition, |^(x)| 2 = 1, we obtain 
the eigenfunctions 



i/>n(x) = \ j\ Sin ( 



wkx 



and the corresponding eigenvalues 

_ h 2 w 2 n 2 
E„ 



2L 2 m 



n = 1, 2, 3, 



n = 1,2,3,... 



(2.21) 



(2.22) 



2.3.2 Harmonic Oscillator 

We can now consider a particle in a classic spring like potential, 

V(x) = -kx 2 . 



(2.23) 
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Figure 2.2: Particle in a box wavefunctions. 



The time-independent Schrodinger equation with this potential is 

h 2 d 2 i){x) 1 2 



2m dx 2 + 2** iKx) = EiKx), 



(2.24) 



we can note that if a reduced mass a = mi . m2 is used we can model the 

' mi+ra2 

behaviour of a chemical bond between two atoms of mass m\ and m<i- A 
simple solution to this Schrodinger equation is given by the fact that as the 
derivative of the wavefunction must give back the square of x plus a constant 
times the original function, the solution takes the form 



^(x) = Ce~ ax2/2 



(2.25) 



However, the most general normalized form of the solution is 

n = 1,2,3, ... 



Mx) = - /4 ^L=H n (y)e-y 2 / 2 



7T \/2 n n\ 
with the energy eigenvalues 



(2.26) 



E n = ftw(n + 1/2), 



(2.27) 



where y = \fax } a — and H n (y) are the Hermite polynomials given in 
Fig. 2.3. These quantum harmonic oscillator states are shown in Fig. 2.4. 
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Figure 2.3: Hermite polynomials and the corresponding energy eigenstates. 




Figure 2.4: Quantum harmonic oscillator eigenstates. 
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2.4 Finite-Difference Discretization 
2.4.1 Time Independent Problems 

In the case for complicated potential fields, V(x) and particle scattering 
models the numerical finite difference method has been used for many years 
to solve the Schrodinger equation [4]. For the time-independent case we can 
simply discretise the Schrodinger equation and put it into matrix form, which 
can then be numerically solved. For the one dimensional case, and ignoring 
the potential, the Schrodinger equation at each point along x can be written 
as 



Er/> Xn = 



h 2 (d 2 ^ 



2m V dx 2 



(2.28) 



Xn 



Now using the basic finite-difference approximation, 

d 2 jj(x) = ip Xn+1 - 2t[) Xn + ^ Xn _ x 
dx 2 x=x n a 2 

where a is the spatial interval spacing, we can write Eqn. (2.28) as 

Eij) Xn = k (2^j Xn - ^ Xn+1 - il> Xn -i) , 



(2.29) 



(2.30) 



where k = 



h 2 

2ma 2 



. This can now be written in matrix form as 





/ 


ipi 


\ 




( 2k 


-k 











... \ 




( 


Ipl 


\ 












-k 


2k 


-k 


















E 




V>3 









-k 


2k 


—k 

























































-k 


2k 


-k 
















) 




v ••• 











-k 


2k j 




V 




J 



(2.31) 



This can also be written in operator form as 



(2.32) 
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where I is the identity matrix. This eigenvalue problem can be solved numer- 
ically, and the corresponding eigenvectors, which represent the eigenstate of 
the particle, can be found. In order to implement a potential, V{x), we can 
simply add a diagonal matrix V to H, where the diagonal components of V 
are equal to the potential at the nodes: V nn = V{x n ) 2 . 

2.4.2 Time-Dependent Problems 
Explicit Method 

The finite-difference discretization of the time-dependent Schrodinger equa- 
tion can be simply done using the explicit method. As before we can discretise 
the spatial part of Eqn. (2.5) using the approximation in Eqn. (2.29). Then 
applying the explicit time-difference approximation, 



dt 



= ^ - ^\ (2.33) 



where b is the temporal interval spacing, we are able to construct the explicit 
finite-difference approximation to the Schrodinger equation: 

h 2 



t t ib 

' X n ' Xn 



2ma 2 Xn+1 n ' rx 'i-i/ ' ' !; x n 

In operator form this can be written 



(2.34) 



= (l-^bH^j ^ , (2.35) 

where as before H is the discretized Hamiltonian (with the potential matrix 
V absorbed) and I is the unit matrix. The problem with this approach is that 
it is numerically unstable and also, more importantly, the operator I — |6H 



2 In finite element analysis, rather than taking the nodal values of the potential, the 
average over the element is taken. 
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is not unitary, which is a required property in order to conserve probability, 
i.e. f ip*ipdx 1. 

Implicit Methods 

Now conducting an implicit discretization we have 



ib 

li 



h 2 



2ma 2 



(2.36) 



Which can also be put into operator form as 



I + i 6H 



(2.37) 



Even though this numerical solution is stable it still does not correspond to 
a unitary transformation, and thus leads to unphysical quantum results. 

Cayley's Form 

A numerical finite-difference technique that produces a stable and unitary 
discretized operator is called the Cayley's Form. For this we use a centered- 
time-difference or Crank-Nicolson Scheme to construct the temporal dis- 
cretization: 



in 



+ 



h 2 d 2 ^(x,t) 



2m dx 2 
h 2 d 2 ^j{x ) t) 



+ V(x)^(x,t) 



2m dx 2 
After spatial discretization we have 

h 2 



+ V(x)4>(x,t) 



(2.38) 



ih 



- 4l 
b 



+ 



2ma 2 

h 2 



2ma 2 



(2.39) 
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After rearranging we obtain 



(2.40) 



where /,„ = ^V^„ and p = Simplifying further we have 

(1 + if Xn + 2ig)^ - ig^ - ig^ 
= (1 - if Xn - 2ig)^J n + ig^y + ig^i n _, , 



(2.41) 



This can then be put in matrix form as 

/ ^ Y 3+1 

(I + «H) ^ 

\ ^ 7 



I - iH) 



As before I is the unit matrix, but now H is given by 
/ A + 2g -g 

-g f2 + 2g -g 

-9 h + 2gk -g 



H = 



V 








-9 /jv-i + 2g -g 
-g f N + 2g ) 



So we have the numerical difference equation in the Cayley's form: 

I-iH 



(2.42) 



(2.43) 



(2.44) 



The temporal operator that relates ip^ to ^■ 7+1 is now not only numerically 
stable but also unitary; this can simply be shown as 

I - iU\* /I-iH\ n + /I-iH 



I + iH/ \I + iH 



I-iH I \I + iH 



= I. 



(2.45) 
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Through this unitary property Eqn. (2.44) then satisfies conservation of prob- 
ability as required, 

^2 ^ +1 *^ +1 a = a = 1. (2.46) 

where a is, as before, the spatial interval spacing Ax. 3 



3 As the spatial axis is discretized we are using summation instead of integration, how- 
ever this is equivalent to the continuum equation given in Eqn. (2.1). 



Chapter 3 

Finite Element Analysis 



Development of the Finite Element Method (FEM) can be traced back to the 
1940's. However, it wasn't until the late 1950's and 1960's that it emerged as 
a useful tool in engineering. Then, when a rigorous mathematical foundation 
was developed in the early 1970's it became a dominant method in applied 
mathematics for numerical modelling of physical systems in many engineering 
and scientific disciplines, e.g. electromagnetic and fluid dynamics as well as 
civil and aeronautical engineering [5]. Olek Zienkiewicz, from University of 
Wales Swansea, originally an expert in finite difference methods (FDM) was 
one of the pioneers in bringing FEM to the wider scientific and engineering 
community through the first book on the subject [6]. 

Even though FEM is a little more complicated to implement compared 
to FDM, one of its biggest advantages is its ability to handle complicated ge- 
ometries (and boundaries) with relative ease. However, even though handling 
complex geometries in FEM is theoretically straight forward, the problem of 
computational time is strongly influenced by the ability to precondition the 
problem i.e. by choosing the most appropriate element type for the most 
efficient computational performance. 

20 
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As a sideline we can also note that FDM is a subset of the FEM ap- 
proach. This can be seen through choosing basis (shape) functions as either 
piecewise constant or Dirac delta function; then the stiffness matrix K can 
be interpreted as a difference operator [7]. Then, by using a uniform mesh 
the FE equations reduce to FD equations. 

There are two specific techniques for the application of FEM to a prob- 
lem, the variational and Galerkin. The variational approach requires a FE 
discretization of the functional associated with the problem (or, if it can be 
defined, the Lagrangian of a system). The discretization is done in the stan- 
dard way using basis functions for each element of the domain considered. 
Then, by minimising the discretized functional and assembling the system 
for all the individual elements we are able to obtain the required FE equation 
of the system. This is a powerful method as it takes into consideration the 
physics of a system 1 in order to simplify and solve the problem. 

As opposed to the variational method the Galerkin approach is directly 
applied to the differential equations of the problem, and then the equation 
can be discretized and assembled in order to obtain the FE equations of the 
system. In this work we will use the Galerkin method as this eliminates the 
work of finding the functional associated with the problem. 

In the rest of this chapter we will lay the foundation for the application 

of FEM to the types of problems we will encounter in the quantum context. 

Beginning with a simple example of a one-dimensional eigenvalue problem, 

we will then go on to discuss time-dependent problems, i.e. the use of a 

combination of FEM, to solve for the spatial part, and FDM, to solve for 
1 Lagrangian of a system for conservative systems or virtual work for the general case 
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the temporal part. Finally, we will give a brief introduction to the use of 
space-time FEM to solve spatial and temporal parts of a problem together. 

3.1 Eigenvalue Problems 

Even though in this work we will be dealing solely with time-dependent prob- 
lems we will never the less include a brief summary of the general eigenvalue 
problem and the use of FEM for their solution. In quantum mechanics the 
eigenvalue problem is one of the most important aspects: in order to de- 
termine energy levels and the associated eigen-functions. For the case of 
complicated potentials, as in irregular lattices, and in quantum dots, numer- 
ical computation of eigenfunctions and eigenvalues is of great importance. 

3.1.1 FEM for Eigenvalue Problems 

To lay down the method of the Galerkin approach to eigenvalue FE problems 
we will consider the case of torsional vibrations of a uniform circular-cross- 
section [8]. The differential equations and boundary conditions required to 
determine the mode shapes and natural frequencies are 












0. 




Rearranging this we can put it into an eigen- value form 



+ 7^ = 



(3.2) 



dx 2 
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2 

where 7 = To apply the Galerkin method we multiply Eqn. (3.2) by a 
test function (j) and integrate it by parts, 



J ^ + = 



— \— —cj)\dx— [ -^--^-dx + [ ^yipcpdx 



q dx I dx J J q dx dx 

d± L _ [ L dAdJ, dx+ f = 

dx y dx y 

=0 

Therefore, once we eliminate the first term on the LHS (due to boundary 
conditions) we obtain 

- / ^^dx+ [ ~fip(j)dx = (3.4) 
J dx dx J 

The next step is to implement a FE approximation using a set of basis func- 
tions, Nf. 

i; = J2ANi = Nij <P = J2 <t>i^i = (3.5) 

i i 

where the basis operator, N, in one dimension is given as 

N = [N 1 ,N 2 ] ^1 = ^ N ^ =] ~r L ( 3 - 6 ) 
The coordinate transformation is then simply given as 

x = ^x i N i . ( 3 - 7 ) 
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Using this information we can write the derivative transformation as 2 

dx X2 — Xi l P 



dtp 

~dj 

d 

dx 



2 

dtp dx 
dx dl; 
2 d 
T e d£' 



2 



dtp 2 dtp 



/ g dec 



(3.1 



(3.9) 



Or in matrix notation we can write 



dtp 



"22 



^2 



(3.10) 



d_ _ I 

dx I p 



Combining this with Eqn. (3.9) we have 

-1 1]=B. (3.11) 

This B operator then gives the FE approximation of the derivative of a 
function ip: 

g = i[-l l]f = B* (3.12) 

where is the FE approximation vector of the continuous function ip. Now, 
using Eqns. (3.12) and (3.5) in (3.4) we obtain 

-E/^^^I^+Et//^^ = (3 - 13) 





' -1 " 




1 



--E 



[-1 1« + 

[Ni N 2 ]tpdi = 



(3.14) 



2 In this ID case the derivative transformation is very basic, however when we go on 
to work in 2D (for space-time FEM) we will require the more complicated 2D coordinate 
Jacobian. 
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where the sum is taken over all the elements. The first term can be simplified 



as 



E-f** 



1 -1 
-1 1 



(3.15) 



and the second term, after a little more computation, can be written as 

^2 J_i [ N 2 Nt N 2 N 2 J y ^ 



e 

The total FE model is then 



2 1 
1 2 



(3.16) 



{ 


1 -1 " 


l 2 7 


'21" 


} 




-1 1 


6 


1 2 





^ = 



(3.17) 



Now, as this will be true for all test functions 0, we can write the two element 
approximation as 



1 -1 




" 2 


1 


" 






-1 2 -1 


- A 


1 


4 


1 


|^ = 


(3.18) 


-1 1 







1 


2 







where A = ^ is the eigenvalue, and using boundary conditions tp = [0 tpi fo] 
is the eigenvector. This matrix torsional vibration equation, with two ele- 
ments, is now in the form of a generalized eigenvalue problem. The com- 
plexities in solving this equation will be dealt with when dealing with such 
problems in the context of quantum mechanics. Just to note, if we wanted 
to solve Eqn. (3.18) analytically we would write it as 



2 - 4A -1 - A 
-1-A 1-2A 



^2 



(3.19) 
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where we have eliminated the first row and column as ip = 0. Therefore, a 
solution exists when the determinant of H vanishes. In this way we obtain 
the eigenvalues and eigenvectors of the problem. 3 

3.1.2 Application to the Schrodinger Equation 

In quantum mechanics the very basic eigen- value problem consists of solving 

h 2 d 2 ^ 



2m dx 2 



V(x)ip = E^j. (3.20) 



The aim is to determine the energy level configurations for particles in various 
potentials and spaces. We will first consider a particle in an infinite well where 
V = 0, this will then be extended to a general well V{x). 

Model of Infinite Potential Well 

To model an infinite potential well using FEM we begin with the Schrodinger 
equation with zero potential, 

^+1* = (3.21) 

where 7 = Now using the FEM construction described in Sec. 3.1.1 we 
obtain the Schrodinger equation FE approximation: 



E 





1 -1 " 


ih 


"21" 




{ 


-1 1 


6 


1 2 


) 



ip = (3.22) 



From the boundary conditions of an infinite potential well we know that the 
nodal approximations at the edges of the well are zero: ipi = ^+1 — 
(where TV is the number of elements). Therefore, assembling for TV elements 



3 For a full solution and explanation of this method see [8] 
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" 2 -1 


••• 


" 




" 4 


1 





... o " 






lp2 


< 


-1 2 
-1 


-1 

2 





ih 

6 


1 




4 
1 


1 

4 







> 


^4 








-1 








'•• 1 










o ••• 


-1 


2 




. 







1 4 _ 


> 







= 



(3.23) 



Model of a General Potential V(x) 

To take into account a general potential we need to model 

d 2 ^ 



dx 2 



v{x)i) + = o. 



(3.24) 



The extra potential term is incorporated through the following FE approxi- 
mation 



2 1 
1 2 



ft 2 Vj-i V ^ & 

e x 7 e 

where V e is the average potential within the element e. Assembling for TV 
elements we obtain the full Schrodinger FE approximation 4 : 

mil l2 ' 







A'+^B'-^C'^ 



(3.26) 



4 An in-depth study of the formulation and solution of quantum eigen- value problems 
can be found in [12]. 
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where the matrices are given as: 



A' 



B' 



a = 



1 -1 

-1 2 -1 
0-12 








-1 1 



V -Vi 







o -v 2 v 2 + v 3 



2 1 





'•• -v N 

-v N V N 



1 4 1 
1 4 



1 2 



and the nodal vector is 



1pN+l 



(3.27) 



(3.28) 



3.2 Time- Dependent Problems 



Before discussing space-time FEM we will first give a basic example of the so- 
lution of time-dependent problems using FEM/FDM. In this example we will 
use the one-dimensional diffusion equation with initial boundary conditions, 
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which take the form 

d 2 ib dib 

kA-^ = pc P A-^- 0<x<L,0<t 

ip(0,t) = % 0<t 

ip(L,t) = 0<t 

ip{x,0) = 0<x<L 



(3.29) 



So the problem we will solve can be written as 



where 7 = We will solve the spatial part of this problem using the 

previous FEM approach, but then the temporal part will be dealt with using 
FDM approach. 

3.2.1 FEM Spatial Discretization 

As in the case of torsional vibrations we begin by multiplying Eqn. (3.30) 
with a test function (j) and then integrating by parts 



70— 
dx 



o 



Jo dxdx J dt 



- 1 So fxfx dX = So ^ (3 - 31) 



Now, using Eqns. (3.12) and (3.5) in (3.31) we obtain 

/ ^ f 1 -J 



- 7 j \ ^B^^de = E S_ x ^ tNtN ^ ( 3 - 32 ) 
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where the sum is again over all the elements, and the vector ip is now time- 
dependent. Simplifying and integrating we obtain 



-72- 2> 



1 -1 
-1 1 



2/3 1/3 
1/3 2/3 



^ (3.33) 



Therefore, we have 

E*M A 



1 -1 
-1 1 



2 1 
1 2 



=0, 



(3.34) 



where A = |f . Then, as this is true for all test functions, 0, we have the 
element equation 



A 



1 -1 
-1 1 



2 1 
1 2 



-0 = 



(3.35) 



When assembled for four elements, we obtain 



A 



1-10 
-12-10 
0-12-10 
0-12-1 
0-12 



2 10 
14 10 
14 10 
14 1 
1 4 



i> = 



(3.36) 



From the initial conditions, in Eqns. (3.29), we know that the first component 
of ^ is a constant (ipo = constant) and the last component is always zero 
= 0). Using this information we can reduce Eqn. (3.36) to 





2 -1 







' 4 


1 


" 




A-00 


A 


-1 2 


-1 


i) + 


1 


4 


1 


j> = 







-1 


2 







1 


4 








where ip = [ipi ip2 ^3] • Writing this in a more convenient notation we have 

AA-i/5 + Bip = C (3.37) 
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3.2.2 FD Time-Integration 

The next step is to use FD techniques to carry out time-integration. This 
can be done in many ways, but here we will only consider two techniques: 
the explicit Euler and the implicit Crank-Nicolson methods. 

Explicit Euler Time-Integration 

This method is the simplest to implement, however it can be unstable. We 
begin with the following approximation 

i, = - (3.38) 

Implementing this into Eqn. (3.37) and simplifying we have 

AtXAip 71 + B [ijj n+1 - i> n ] = AtC 
[A / A / - B'] ^ n + B f ip n+1 = C (3.39) 

where A / = A, B / = B, C' = AiC and A / = AiA. Using the initial condition 
vector ijj we can determine = B'^/r, which can be solved to obtain ip 1 . 
This process can be continued in order to obtain ip 2 from ip 1 and so on. 

Implicit Crank-Nicolson Time-Integration 

The explicit Euler method is mathematically and computationally very sim- 
ple, however it can be unstable. On the other hand, the implicit Crank- 
Nicolson method in unconditionally stable, even though it is slightly more 
complicated and computationally intensive. In order to use this method we 
begin with the following approximation 



~ _ 1 f r ~ ~ _"| n r „ _ _i n+1 ^ 

B^ = - [C - XA^ + [ c - AA ^J \ (3-40) 
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+1 3 




-1 
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+1 

-1 2 





X 



Ax 



x 



Figure 3.1: Node numbering and coordinates of rectangular element. 



Therefore we have 



B^ n+1 - B^ n 
[B 7 + A 7 A 7 ] 



At + ~ n + r 

c 'n , c 'n+l 



2 2 
[B 7 - A 7 A 7 ] ^ n 



(3.41) 



where A 7 = A, B 7 = B, C 7 = and A 7 = ^ . As stated previously 
this method is unconditionally stable; so even though oscillations occur and 
the accuracy may suffer for large step sizes At, the oscillations never become 
unbounded. 5 . 

3.3 Space-Time Finite Element Method 

In order to solve Eqn. (3.30) using space-time FEM the x, t domain will be 
discretised into rectangular elements, labeled as in Fig. 3.1. The approximate 
solution can then be modelled using the linear, dimensionless, local basis 



3 For further details and example of time-integration techniques see [8] 
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functions: 



(1-x)(1-t) 



N 2 (x,t) 



(1 + x)(1-t) 



4 



4 



W 3 (x,r) 



(l + x)(l + r) 



a^(x,t) 



(l-x)(l + r) 



(3.42) 



4 



4 



The other important factor to consider is that of continuous or discontinuous 
boundaries between the temporal elements t n and t n+ i. In the discontinuous 
method the approximate solution of ip is continuous within the elements but 
discontinuous at the boundaries, so each time step forms a slab in space. In 
the continuous case the solution flows uninterrupted from one time step to 
another, like a continuous surface in the x, t domain. It has been shown (in 
[10]) that the linear discontinuous method has a higher accuracy compared 
to the continuous method. This is due to the extra nodal degree of freedom 
available in each time step. However, in order to implement the discon- 
tinuous method an extra jump term has to be included in the discretization 
process. To demonstrate the space-time FEM process we will show the linear 
continuous discretization, and then demonstrate the discontinuous case. 

3.3.1 Linear Continuous Discretization 

We begin by writing the solution as a function of the local coordinates 



4 




(3.43) 



i 
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where TV [TVi, 7V 2 , Af 3 , 7V 4 ] and ip [^1,^2,^3,^4]. The corresponding 
differentials are 



dx 
dip 

Writing as a matrix this becomes 



dip dx dip dt 

dx dx dt dx 

dip dx dip dt 

dx dr dt dr 



(3.44) 



- dip - 




dx 


dt ' 




" d± ' 




■ dip ~ 


% 




ox 






dip 


= J 


% 


L dr J 




L dr 


dr J 




L dt J 




- dt - 



(3.45) 

where J is the Jacobian of the transformation. Using the coordinate trans- 
formations 

4 4 

x = Y^ x i N i = Nx t = Y^ f i N i = Nt (3.46) 
1 1 

and the nodal numbering described in Fig. 3.1, the Jacobian reduces to the 
simple form 



J = 



— 

2 ; 

At 







(3.47) 



The inverse of this, which will be required later, is simply given as 



j 1 



h 
h 



(3.48) 



Therefore we now have 



- dip - 




- dip ~ 




= J 1 


% 


L dt J 




L dr J 



(3.49) 



The local derivatives with respect to x an d t can be discretised as 
dip 1 



dx 



[^(1 - r) - ^ 2 (1 - r) - ^(1 + r) + ^{1 + r)] 



^ = -^[^i(l-x)+^(l + x)-^3(l + x)-^4(l-x)] (3-50) 
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In matrix form this can be written as 



- dip - 


1 




% 

L dr J 


~ ~4 



(l_ r ) _(l_ r ) _(l +r ) (l +r ) 

(i-x) (i + x) -(i + x) -(i-x) 



^3 
^4 



Now, combining Eqns. (3.49) and (3.51) we have 

d± 



■ Aip. 
(3.51) 

(3.52) 

(3.53) 

At At ~ 

In order to use this in the discretisation of the diffusion equation we write 



= J~ l Aip = Bip, 



where 



B = 



d± 

at 



(i-t) 

Ax 

(i-x) 
At 



Ax 
(1+X) 
At 



-(l+r) 
Ax 

-(i+x) 



(l±l) 
Ax 

-(i-x) 



(3.54) 



where Bx and B 2 are the upper and lower rows of B. The continuous space- 
time discretisation begins with Eqn. (3.31), and then we apply the discrete 
derivative and function approximations to obtain 6 

r+1 ' - AxAt 



/it 



£ (NtB 2 + 7 B{B 1 )^dxdr 



1> 



(3.55) 



Carrying out the simple integrals, and noting that this is true for all test 
functions, we obtain the following set of equations 7 



= 



2 1-1-2 
12-2-1 
12-2-1 
2 1-1-2 



2jAt 
Ax 2 



2-2-11 


1 






-2 2 1-1 






^2 


-11 2-2 




1 




1-1-2 2 






. ^4 



(3.56) 



6 Note that as we are transforming from a global (x, y) coordinate basis to a local (%, r) 
basis we have implemented the volume transformation dxdt = \ J\dxdr. 

7 Where the extra factor of 2 in the second term comes from integrating over r i.e. 
J-i dr = 2, even though there are no explicit r variables. 
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Taking advantage of boundary conditions we can note that for each time step 
t n to £ n+ i the nodal values for t n are known. In this way we can reduce the 
4x4 element equation to a 2 x 2 one. However, in order to pick the correct 
set of equations in Eqn. (3.56) we need to look at the structure of the test 
function: 



Ncj) = [ Ni N 2 N 3 N 4 ] 



(3.57) 



Here we can note that the shape functions for Ni and N 2 are 1 at time i n , and 
N$ and N± are 1 at time t n +\. Thus if we use the first two rows in Eqn. (3.56) 
(associated with N\ and 7V 2 ) we obtain an explicit numerical method which 
is weighted on information from t n , however if we use the second two rows 
(associated with A3 and N±) we obtain an implicit numerical method which 
is weighted on information from t n +\. Thus, going for the second two rows, 
and rearranging the results we have 

= £ 





"12" 




" -1 1 


) 


^3 


( 


2 1 


+ Ax 2 


1 -1 




. ^4 _ 



! 


"21" 


27At 


1 -1 " 


1 






1 2 


Ax 2 


-1 1 




. ^2 _ 



(3.58) 

A space-time difference stencil can now be obtained by assembling the above 
equations for two neighbouring elements, which can then be easily extended 
for larger number of elements. 8 



B For a detailed discussion and analysis of space-time FEM see [10] 
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Figure 3.2: Linear discontinuous finite elements in space-time. 

3.3.2 Linear Discontinuous Discretization 

For this method the solution is linear within each time step and discontinuous 
at the temporal boundaries, Fig. 3.2. In order to implement this method we 
need to include a jump term, 



into the space-time discretization process. To begin the discontinuous space- 
time discretization we first write the approximate element solution for the 
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time interval t n to t n+1 as 

ij) = Nip = [ Nt N 4 N 2 N 3 ] 



^3 



— [ Nj >n N jtTl+1 N j+lj1l N j+1<n+1 ] 



, (3.60) 



^j+l,n+l 

where the shape functions N are as in Eqn. (3.42), and the relabelling is 
done to make the final assembly process simpler (to take account of this 
relabelling we will also need to rearrange the rows and columns of Eqn. (3.56) 
for the discontinuous space-time element equation). Using this relabelled 
approximate solution we can discretise the jump term as 



/ </> +t iV +t [N + i> + - N'ij-) dx 

( 





N- 

















[Nj, n N j+1 , n 0] 



V 



[0 Nj, n N j+1 , n ] 



0+ 

^j+l,n+l 









fin 
V'j+l.n-l 








) 





(3.61) 



Calculating the tensor products, and for t n we set r = —1 in AT, we have 

4> 
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o (i + 0(i + 6 o 
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(3.62) 
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Carrying out the integrals we obtain 
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^j+l,n-l 



(3.63) 



We now multiply this by — ^ so it can be added to Eqn. (3.56). However, we 
must also rearrange the matrices in Eqn. (3.56) so that they correspond to the 
nodal ordering in Eqn. (3.60). After doing this we have the full discontinuous 
space-time element equation 
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. rj+l,n+l 



















fin 

(3.64) 



Writing this in operator form we have 



(A' + B' + C) 











= D' 


V'i.n 




V'j+l.n-l 


. ^'+l,n+l . 







(3.65) 



where A', B', C', and D' are the operators of the respective matrices in 
Eqn. (3.64). 



Chapter 4 

Time-Dependent Analysis 



In this chapter we will begin by discretising the time-dependent Schrodinger 
equation by the use of FE for the spatial part and Crank-Nicolson for the 
temporal part. As there are extensive results and literature on this method 
([11] and references therein) we will have a comparable benchmark. We will 
first do this analysis for a Gaussian wave-packet in an infinite potential well, 
and then we will conduct a similar analysis but with a potential barrier lo- 
cated within the well. The next step will be to discretise the time-dependent 
Schrodinger equation using the space-time FE approximation. The results of 
the space-time method can then be compared to those of the first method. 

4.1 Crank-Nicolson and Finite Element Analy- 
sis 

4.1.1 Infinite Potential Well 

We will begin with the simple case of a particle in an infinite well, Fig. 2.1. 
This will then form the basis for modelling a particle in an infinite well with 
a finite potential barrier. 



40 



CHAPTER 4. TIME-DEPENDENT ANALYSIS 



Crank-Nicolson Temporal Approximation 



We begin by discretising the equation 

dip .^d 2 ip 
dt dx 2 

where b — Now, following the steps of Sec. 3.2, we first apply 
discretisation: 



r +1 

J N f Nd^ 



h r +1 

+ i — / B f Bd^ = 

2m J_ x 



After computing the integrals we obtain 



2 1 
1 2 



1 -1 
-1 1 



V» = o. 



In operator form we have 



A-0 + zB^ = 0, 



where 



A = 



2 1 
1 2 



and B = 



6k 
2mll 



1 -1 
-1 1 



Now applying the Crank-Nicolson approximation we obtain 



A^ n+1 - A^ n = - ( — iB^ n + — iB^ n+l 

2 2 



and rearranging we have 



(A + iB) V n+1 = (A - iB) ip n . 



The constants ^ have been absorbed into B, which gives 



B 



Aml 2 e 



1 -1 
-1 1 
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Here At is temporal difference and l e is the spatial element size, m is the 
mass of the particle and h is the Planck constant. It can also be seen that 
the transformation in Eqn. (4.7) is unitary, and as an extra confirmation of 
its validity it takes the same form as the full finite difference approximation 
in Eqn. (2.42). 

Construction of the Numerical Method 

Our aim now is to solve Eqn. (4.7) for ip n+1 given ip n i.e. find ip(t) knowing 
the initial condition ^(0). In order to conduct this iterative computational 
calculation we must first simplify Eqn. (4.7) so that the complex values can 
be easily handled. The element state vector ijj currently takes the form 



'A ' 




" a 1 + ib 1 ' 






a 2 + ib 2 



(4.9) 



where ipi is the complex left nodal value and ^ 2 is the complex right nodal 
value. However, we can write this complex two-component vector as a real 
four-component vector: 







' a 1 ' 


Im[ipi] 




b 1 






a 2 


_ Im[ip 2 ] 




b 2 



(4.10) 



Using this real four-component element vector we can write the complex 
element equations in a totally real form as 



(A' + aB') ij n+1 = (A' - qB') 4> n , 



(4.11) 



where a = §^§1, and the real matrices are given as 



A' 



2 10 
2 1 
10 2 
10 2 



B' 



0-101 
10-10 
10-1 
-10 10 



(4.12) 
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Eqn. (4.11) can now be assembled using the normal FE method. For example, 
considering two elements the nodal vector becomes 



ij T = [ a 1 b 1 



a 



2 Ir a 3 b 3 ], 



(4.13) 



and the matrices take the form: 



A' = 



2 1 

1 

1 
10 



2 






4 



1 

1 

4 1 

2 

1 2 



B' = 




1 


-1 








1 









-1 


2 


-1 



1 


-2 

1 








-1 



1 






1 





(4.14) 



Initial State Function 



At the initial time t = we can assume that a particle is placed into an 
infinite potential well at position x with a momentum fc . The initial par- 
ticle state can then be modelled as a Gaussian wave packet, as described in 
Appendix A. 2: 



ip(x) 



Re[ip(x)] 
Im[ip(x)] 



(4.15) 



The value of this wave packet at each nodal position j can then be written 

as 



Re[^} 
Im[i/jj] 



(j- 

V27R7 2 



-(xj-x ) 2 /4a 2 



cos(/c Xj) 
sin(/c Xj) 



(4.16) 



Doing this for each node of the spatial domain we can construct the inital 
state vector: 



(^°) T = [ ifc#x] Re[ip 2 ] Im[tp 2 ] 



Re[ip n+1 ] Im[ip n+1 ] ] , 
(4.17) 
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where n is the number of elements. 

We can note that if we use a = 2 the wave-packet in Eqn. (4.16) is 
automatically normalized: 



This then removes the added task of normalizing the final results. 
Numerical Solution 

We can now determine the time evolution of the state vector in Eqn. (4.17). 
We begin by first finding 



Once we evaluate this simple matrix multiplication we obtain the vector 



which can be achieved by using the simple LU decomposition method. Once 
the solution for ip 1 is obtained, we repeat the process to find and in turn 
i/j 2 and so on until we reach the solution for ip at time t. 

Numerical Results 

After running the numerical simulations (code described in Appendix B) 
with the parameters: 250 elements, time-step dt = 0.5, infinite well size of 
—20 < x < 20, and the wave packet initially centered at x = 0, we obtain 
data for the real and imaginary parts of the time evolution of the wave packet. 




(4.18) 



4>° = (A' - oiB') tp°. 



(4.19) 
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A selection of these results and their square-sums (Re 2 + Im 2 ) are plotted in 
Figs 4.1 and 4.2. In these plots it can be seen that the wave packet moves 
to the right until it collides with the infinite potential barrier on the right. 
It is reflected, and then it continues to the left side of the well, where it 
again rebounds to head back to the right. On collision with the infinite walls 
the Gaussian envelope undergoes a distortion. This is due to the fact that 
the real and imaginary parts of the wave-packet, even though they are not 
physically observable, undergo phase changes on reflection. If this simulation 
is run long enough the Gaussian envelope will spread out until it covers the 
entire well [2] 1 . 

In order to study the conservation of probability we incorporated a simple 
trapezoidal rule to sum the areas under the elements at each time step - this 
was done for varying element sizes in order to study the accuracy. In Fig. 4.3 
it can be seen that when 50 elements or more were used the total area after 
each time step remains constant at 1, as expected. Using only 5 elements 
the area drops to « 0.43193, however it remains constant at each time step. 
The difference in areas for the varying number of elements can simply be 
accounted for by the fact that with a lower number of elements the shape of 
the solution can not be resolved accurately, hence the area under the curves is 
not representative of the actual area. On the other hand, as the area is always 
constant irrespective of the number of elements, it implies the conservation 

of probability property is maintained. 

1 This can be seen in Fig. 4.12 where a low energy wave-packet is placed in an infinite 
well divided by a finite barrier (which is larger relative to the energy of the packet by a 
factor of 2.5). As the wave-packet is trapped on the left side, and very little is transmitted 
to the right side, it eventually spreads and covers the entire left half of the well. 




Figure 4.1: Wave packet trapped in a infinite well (cont. on next page). 




Figure 4.2: (Cont. from previous page) Wave packet trapped in an infinite 
well. 
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Figure 4.3: Conservation of area at each time step for varying number of 
elements. 



Wave- Packet with k = 

If we set the packet wave-number, fc , in Eqn. (4.16) to zero, we have a 
stationary wave-packet, which represents a particle at rest. The numerical 
results in this case behave as expected, as the packet disperses and spreads 
over the infinite well, Fig. (4.4). We can also note that as the packet spreads 
over the well with time, the area remains constant, Fig. (4.5). 

4.1.2 Infinite Potential Well with Barrier 

To model the infinite potential well with a barrier, Fig. 4.6, we will discretise 
the Schrodinger equation with the potential term: 

dib h d 2 ib i 
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Figure 4.4: Wave packet with k = spreads over the infinite well (Re 2 + Im 2 
are shown). 
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Figure 4.5: Conservation of area at each time step for packet with fc 
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Figure 4.6: Infinite potential well with a barrier. 
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where V is a constant potential for the barrier and zero everywhere else in 
the well. The LHS of Eqn. (4.21) will be discretised as before, so we will 
only consider the element equation for the RHS. Following the procedure of 
Sec. 3.1, due to the similarity of the potential term and the eigenvalue term, 
we derive the potential element equation through the Galerkin technique as 



e 







fj 





[Ni N 2 ]$<% 



2 1 
1 2 



(4.22) 



We can then write the total Schrodinger element equation, including the 
potential term, as 



2 1 
1 2 



7 • f 6h 



1 -1 
-1 1 



2 1 
1 2 



V» = o. 



(4.23) 



Writing this in operator form we have 

+ ? |b + c}^ = 0, 



where 
A = 



2 1 
1 2 



B = 



6h 

2ml 2 e L 



1 -1 
-1 1 



and 



(4.24) 



2 1 
1 2 

(4.25) 



Applying the Crank-Nicolson approximation, as before, we obtain 

(a' + i {b' + c'}) = (a' - « {b' + c'}) r , 



(4.26) 



where A = A, but B and C are given as 



_ 6HAt 
4m/g L 



1 -1 
-1 1 



and 



C = 



2h 



2 1 
1 2 



(4.27) 



CHAPTER 4. TIME-DEPENDENT ANALYSIS 52 
Numerical Construction 

In order to compute Eqn. (4.26) numerically we separate the real and complex 
parts as in Sec. 4.1.1. In this way we obtain the four component element 
vector as in Eqn. (4.10), and the element equation becomes 

(A 7 + aB' + (3C) i; n+1 = (A' - oB' - PC) $ n , (4.28) 

where a = and f3 = and the real matrices A' and B / are as in 

Eqn. (4.12), and C is 

"0-20-1 
r , = 2 1 
0-10-2 
_ 1 2 

Using the initial state function as in Sec. 4.1.1 the numerical solution of 
Eqn. (4.28) now follows the same procedure as in Sec. 4.1.1. The most 
important difference is that for the elements corresponding to the potential 
barrier f3 ^ 0, but for all other elements where j3 = we have the original 
element equation. 

Numerical Results 

To implement the barrier modifications to the initial infinite well we begin, as 
before, with an infinite well of size —20 < x < 20, and then we incorporate a 
potential barrier of height V = 2.5 at the center of the well: —0.8 < x < 0.8. 
Therefore, if we consider a total of 250 elements the barrier is located at 
elements 120 < x e < 130; we can then assemble Eqn. (4.28) with V e — 
from elements 1 to 119, then with V e = 2.5 from 120 to 130, then again with 
V e = from 131 to 250. Then using a time step dt = 0.5 and beginning with 



(4.29) 
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the initial wave-packet at x = — 13 (located to the left of the barrier) we 
can obtain the data for the real and imaginary parts of the time evolution 
of the initial wave packet as before. A selection of these results and their 
square-sums are plotted in Figs 4.7 and 4.8. In these plots it can be seen 
that the wave-packet moves to the right until it hits the barrier. A small 
part of it is transmitted through the barrier and the rest is reflected back. 
The reflected and transmitted parts continue moving until they collide with 
the infinite walls of the well and return back to the barrier. 

When we look at the conservation of probability property we find that 
with the introduction of a finite barrier the conservation of total area now 
fluctuates around the previous value of 1, Fig. 4.9. Comparing Figs. 4.7 
and 4.9 (for 250 elements) we can conclude that the first fluctuation (dip at 
5 < t < 10), in Fig. 4.9, occurs when the wave-packet first interacts with the 
barrier, in Fig. 4. 7. The second fluctuation (peak at 15 < t < 25) occurs when 
the reflected and transmitted packets interact with the well walls. However, 
after the second fluctuation (dip at 28 < t < 35, the transmitted and reflected 
waves recombine, in this case when the packet interacts with the well walls 
no peak fluctuation occurs. The next fluctuation is again a dip when the 
packet interacts with the barrier again. These fluctuations then continue, 
with the dips representing an interaction with the barrier and the peaks 
representing interactions with the well walls (when reflected and transmitted 
waves exist). We can argue that when the problem involves reflected and 
transmitted waves we require greater resolving power in order to determine 
the shape of the total probability distribution ip. We can show this to be 
the case by increasing the number of elements in the simulation and keeping 




Figure 4.7: Wave packet trapped in a infinite well with barrier at elements 
120 < x e < 130 (cont. on next page). 
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Figure 4.9: Conservation of area at each time step for varying number of 
elements. (Also on the same plot is area conservation for 500 elements and 
time step dt=0.1) 



the timestep dt constant. In Fig. 4.9 it can be seen that if we increase the 
number of elements the peak and dip fluctuations begin to diminish - for 
150 elements the fluctuations have a maximum value of 2.2%, but for 1000 
elements the fluctuations fall to a maximum value of 0.2%. 

Also, from Fig. 4.10 we can see that if we decrease the timestep (dt) but 
keep the number of elements constant we can slightly decrease the fluctua- 
tions. However, for time-steps smaller than dt = 0.1 there seems to be no 
change in the fluctuations. Taking these observations into account we con- 
ducted the barrier simulation for 500 elements and a time step of dt = 0.1. 
From Figs. 4.9 and 4.10 we can see that this is slightly more effective than 
just using 500 elements (and time step dt = 0.5) or a time-step of 0.1 (and 
250 elements) alone. 
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Figure 4.10: Conservation of area at each time step for varying size of dt. 
(Also on the same plot is area conservation for 500 elements and time step 
dt=0.1) 



In order to study these fluctuations in detail we conducted further simula- 
tions, but this time by varying the initial wave- vector fco, given in Eqn. (4.16). 
All the previous simulations were performed with k = 2, and the finite po- 
tential barrier was of height V = 2.5. So if we used k « V the low energy 
wave packet would be effectively trapped in the left side of the well - with 
very little being transmitted. On the other hand if we used ko » V , the 
high energy wave-packet would move around the infinite well unhindered by 
the barrier - so very little reflection would occur. As these cases will be 
very similar to our first set of results (of the wave-packet in an infinite well) 
the fluctuations due to the interaction with the barrier should disappear. In 
Fig. 4.11 we see that using a value of fc = 3 the majority of the wave-packet 
is transmitted. And in Fig. 4.12 where we used ko — 1 very little is trans- 
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mitted, so the packet is effectively trapped in the left side. When we plot 
the areas using these wave vectors, and also of simulations with fc = 0.1, 
4, and 12 in Fig. 4. 13 we find that with a very small and very high fc the 
fluctuations do indeed decrease (for the case of fc = 0.1 and 12 they almost 
vanish) . From this we can conclude that when the energy of the wave-packet 
is comparable to the potential barrier, and so any interaction between the two 
becomes significant, we require greater number of elements to take account of 
the finer resolution changes in the wave-packets Re and Im components. So 
the fluctuations in probability conservation are more to do with the spatial 
element discretisation rather than the size of the time step dt. 

4.2 Space-Time Finite Element Analysis 

4.2.1 Linear Continuous Space-Time Analysis 
Infinite Potential Well 



In order to apply space-time analysis we can follow the mathematical de- 
scription for the discretisation of the diffusion equation in Sec. 3.3.1. As the 
Schrodinger equation with V = 0, 



~dt 



7 ^ = 



7 = i 



h 

2m' 



(4.30) 



is similar to Eqn. (3.30), we can simply write the discrete space time form of 
the Schrodinger equation as 







2 1 

1 2 

1 2 

2 1 



-1 
-2 
-2 
-1 



-2 
-1 
-1 
-2 



. HAt 

mAx 2 



2-2-11 

-2 2 1-1 

-11 2-2 

1-1-2 2 



^3 
^4 



(4.31) 
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Figure 4.11: Wave packet with initial wave vector £r = 3 trapped in a infinite 
well with barrier at elements 120 < x e < 130. 
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Figure 4.13: Conservation of area at each time step for varying initial wave 
vector fc . 



where ipi and ip 2 are j and j + 1 nodal values respectively at t n and ^4 and 
^3 are j and j + 1 nodal values respectively at t n +i (Fig. 3.1). As ipi and ip 2 
are known at each time step we can split Eqn. (4.31) into two. As we noted 
in Sec. 3.3.1, using the first two rows we obtain an "explicit" time difference 
method, as the shape functions for Ni and N 2 are weighted at t n . After some 
rearrangement we have 



2 1 
1 2 



+ i- 



. hAt 


1 -1 " 




^4 


mAx 2 


-1 1 


1 


. ^3 . 





"21" 


. 2ftAt 


1 -1 " 




"^1 " 


1 


1 2 


mAi 2 


-1 1 


) 


. ^2 _ 



. (4.32) 



We can do the same for the rows associated with the shape functions N 3 
and A^ 4 , giving an "implicit" time difference method, as this time the shape 
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functions are weighted at t n+1 . From these we obtain 

2hM 



2 1 
1 2 



+ i 



mAx 2 



1 -1 " 


) 


V>4 


-1 1 




. ^3 . 



2 1 " 


. fiAt 


1 -1 " 






1 2 


mAx 2 


-1 1 


) 


. ^2 _ 



(4.33) 



These equations can then be written as 



[A + i2B 



A - i2B] -0 n , 
A - iBl ^ n 



where 



A = 



2 1 
1 2 



B = 



hAt 



1 -1 
-1 1 



(4.34) 



(4.35) 



mAx 2 

We can note immediately that both these time step methods do not form 
unitary transformations, hence the conservation of probability property is 
not maintained 2 . In Fig. 4.14 we have plotted the areas at each time step for 
these system of equations, and we can confirm that the area is not conserved. 
For the explicit case the area explodes, and for the implicit case, even though 
it remains finite, it is heavily damped and so it drops to zero. Also, if we 
combine Eqns. (4.34) and (4.34), by summing them, we obtain 



A + i§B 



n+1 



A-i-B 

2 



(4.36) 



which is numerically identical to the equation derived by the FE and Crank- 
Nicolson methods in Eqn. (4.7). The reason for this is that when we add 
the two forms we are averaging the information from t n and t n+1 and so we 
reproduce the Crank-Nicolson method (in particular for the linear elements.). 



2 For a difference method A0 n+1 = Bip n to be unitary it must satisfy ^ n + 1 t^ n + 1 
^ nt (A- 1 5)t(A- 1 5)^ n = ^ n H n 



CHAPTER 4. TIME-DEPENDENT ANALYSIS 



63 




Figure 4.14: Area after each time-step for explicit scheme (left) and implicit 
scheme (right). Area explodes for the explicit case, and remains finite but 
not constant for the implicit case. 



Infinite Potential Well with Barrier 

Even though we have shown that working with the continuous space-time 
method is actually a disadvantage (due to a lack of conservation of proba- 
bility) we will continue to construct the space-time potential term for the 
Schrodinger equation as this will also be implemented into the discontinuous 
method. 

In order to apply the space-time discretisation to the potential term we 
proceed as follows (assuming element sizes, Ax and At, are constant over 
the space-time): 

j r L At At - 

--J V[x)^dx = -i^—^2V e J (t) ] N ] N^d X dr 

= ^ {-*^^ ^ /7 ^^dxd-r j ^, (4.37) 

where the space-time shape vector components, N = [TVi, 7V 2 , Af 3 , A^], are 
as in Eqn. (3.42). After integrating over \ an d t 5 we obtain the element 
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equation 



— i- 



.AxAtV e 
36^. 



4 2 12 
2 4 2 1 
12 4 2 
2 12 4 



(4.38) 



To equate this to Eqn. (4.31) we need to multiply this by the constant — 
this is due to the fact that in Eqn. (4.44) we divided the entire system by 
Ax and multiplied by —12 in order to simplify the system. After doing this 
the continuous space-time element equation is given by 

-1 1 
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,(4.39) 



Now, reconstructing the explicit and implicit parts as before, we obtain 



A + iB + iC 
A + i2B + i2C 



A — i2B — i2C 



A-iB-iC 



(4.40) 
(4.41) 



where we have 

A- 2 1 
A 12 



B = 



hAt 

mAx 2 



1 -1 
-1 1 



3h 



2 1 
1 2 



(4.42) 



Just from observation we can conclude that these are not consistent with the 
conservation of probability. However, as before, when we add Eqns. (4.40) 
and (4.41) together we obtain 

-0 n+1 = 



~ 3~ 3 ~ 
A + i-B + i-C 
2 2 



A-i^B-i-C 
2 2 



(4.43) 



which is numerically identical to Eqn. (4.26). 
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4.2.2 Linear Discontinuous Space-Time Analysis 
Infinite Potential Well 

By following the procedure described in Sec. 3.3.2 we can write Eqn. (4.31) 
in linear discontinuous form as 
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(4.44) 



Separating this into real and imaginary parts, as before, we obtain 8x8 
matrix equation: 

(A' - B' - C) $ a+1 = -D'^, (4.45) 
where the vectors are given as 



tpa+l = 



Im[^ n ] 

Im bPj,n+l\ 



R e bP]J 

Im bPj,n] 

Re[^ +l ,n-i\ 
RebPj+iJ 



(4.46) 
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and the matrices are given as 
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(4.48) 



(4.49) 



(4.50) 



These can then be assembled by summing the lower-right 4x4 components 
of the first element to the top-left 4x4 components of the second element, 
and so on for n elements. We can note from Eqn. (4.45) that this system of 
equations is not unitary. 



Numerical Results 



Using the well specifications as in the Crank-Nicolson method (—20 < x < 
20) we obtained the wave-packet time evolution results. In Figs. 4.15 and 
4.16 it can be seen that even for the simple case, a wave-packet in an infinite 
potential well, there is a small amount of damping in the conservation of 
probability plots. This damping can be controlled by varying the number of 
elements and the size of the time step. Fig. 4.15 shows that for 100 elements 
there is significant damping for time steps dt > 0.1; but for time steps dt < 
0.05 the damping becomes negligible, however such a small time step comes 
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Figure 4.15: Conservation of probability at each time step for 100 elements 
and varying size of timestep size: dt. 

at the cost of greater computing time. In Fig. 4.16 we have shown that the 
damping can also be controlled, to a lesser extent, by increasing the number 
of elements. Using 50 elements we have an almost 2% loss of probability after 
100 time steps, however with > 100 elements this loss reduces to 1%. For n 
elements we must again consider the computation time; for the discontinuous 
space-time method we have 8x8 element matrices which gives global matrices 
of order 4n + 4. Therefore, doubling the number of elements will increase the 
global matrices by a factor of 4 which would in turn require more computation 
time. 

In Figs. 4.17 and 4.18 we have plotted the total probability at timestep 
t = 100 (for 100 elements) against variations in time-step dt. It can be 
seen that for smaller and smaller time steps the damping almost decreases 
to zero. Also, in Fig. 4. 19 we have plotted the total probability at timestep 
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Figure 4.16: Conservation of probability at each time step for dt = 0.05 and 
varying number of elements. 

t n = 100 (for dt = 0.05) against variations in number of elements. This 
shows that changing the number of elements has very little affect after 100 
elements. This implies, for the simple case of a packet in the infinite well 
without any barrier interactions, the time step size dt has more of an effect 
on the damping than the number of elements used. 

In Fig. 4.20 we have made a comparison of the possible number of ele- 
ments, which do not require large computation times, and varying timesteps 
dt. We can see that using 100 elements with dt = 0.05 gives accepptable 
damping of < 1%, where if we use dt = 0.01 we have almost no damping but 
the computation times are considerably larger. 

When we compare the discontinuous space-time method to the Crank- 
Nicolson and Finite-Element method in Sec. 4.1.1, we can see that using the 
Crank-Nicolson method is not only computationally efficient but it also has 
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Figure 4.17: Conservation of probability at t = 100 against dt for 100 ele- 
ments. 




Figure 4.18: Close-up of conservation of probability at t = 100 against dt for 
100 elements. 
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Figure 4.19: Conservation of probability at t = 100 against number of ele- 
ments for dt = 0.05. 




Figure 4.20: Conservation of probability at each time step for 50 and 100 
elements and timesteps dt — 0.1. dt — 0.05, and dt = 0.01. 
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zero damping i.e. it holds the unitary property of the Schrodinger equation. 
Infinite Potential Well with Barrier 

We can include a potential barrier into the discontinuous space-time method 
by using the continuous space-time potential term, given in Eqn. (4.39): 



' 3h 



4 2 2 1 
2 4 12 
2 14 2 
12 2 4 



(4.51) 



where the matrix has been rearranged to take into account the nodal num- 
bering in the discontinuous method. Writing this in real form we have the 



8x8 matrix 



v 



VeAt 

3h 



(4.52) 



We then have the full discontinuous space-time element equation: 

(A' - B / - C' - V) = -D^ a , (4.53) 

where A', B', C and D' are as before in Eqns. (4.47)-(4.50). 
Numerical Results 

As for the Crank-Nicolson method the barrier of height V e = 2.5 is imple- 
mented to the center of the well: —0.8 < x < 0.8; then the assembly process 
is carried out as before: for 250 elements we assemble Eqn. (4.28) with V e — 
from elements 1 to 119, then with V e = 2.5 from 120 to 130, then again with 
V e = from 131 to 250. Then using a timestep dt and beginning with the 
initial wave-packet at x = — 13 we obtain the time evolution of the initial 
wave-packet by the use of LU decomposition. 
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Figure 4.21: Conservation of probability at each timestep, for various number 
of elements (timestep size dt = 0.01). 



In Figs. 4.21 and 4.22 we have the conservation of probability plots for 
variations in number of elements and timestep size. Fig. 4.21 shows that 
by using larger number of elements we can resolve the wavefunction in more 
detail and thus reduce the fluctuations. However, one point of interest is 
the fact that there are no peaks in the fluctuations, as was the case for the 
Crank-Nicolson method, Fig 4.9. 

The main issue with the discontinuous space-time method is the damp- 
ing and the computation time. In Fig. 4.21 the damping can be seen to 
decrease by decreasing the timestep size (dt), however this in turn increases 
the computation time significantly. 

In Fig. 4.23 it can be seen that for the same parameters (250 elements 
and dt = 0.05) the Crank-Nicolson method conserves probability perfectly, 
whereas the discontinuous space-time method suffers from damping and thus 
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Figure 4.22: Conservation of probability at each timestep, for various 
timestep sizes dt (100 elements). 



falls by 1% after 100 timesteps. However, the fluctuations assciated with the 
wavepackets interaction with the barrier and the well walls are slightly less in 
the space-time method compared to the Crank-Nicolson method. In terms of 
computation time the Crank-Nicolson method is almost 5 times faster than 
the space-time method for equivalent parameters. 
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Figure 4.23: Conservation of probability at each timestep, for discontinuous 
space-time (DST) and Crank-Nicolson (CN) methods for 250 elements and 
timestep dt = 0.05. 



Chapter 5 



Wave-Packet in Sinusoidal 
Potential 



The Crank-Nicolson method turned out to be more efficient and accurate 
compared to the space-time methods. An important limitation of the Crank- 
Nicolson time evolution method was when we model the wave function in- 
teracting with a change of potential. This interaction results in fluctuations 
in the conservation of probability, however it was seen that these fluctua- 
tions can be controlled by increasing the number of spatial elements used. In 
this chapter we will attempt to test this limitation by applying the Crank- 
Nicolson method to a particle in a periodic lattice potential, 



where A is the potential amplitude and k is the wave number. This periodic 

lattice potential is commonly used in solid state physics in order to simulate 

the atomic lattice structure within materials 1 , Fig. 5.1. In terms of boundary 

conditions we could use periodic conditions, ip(x = 0) = i/j(x = L), however 

in order to implement the changes into our previous work we will continue 
1 For a detailed study of the solid state application and simplified solutions see [13] 



V(x) = A cos [kx] , 
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Figure 5.1: Model of a periodic lattice potential within a quantum wire. 

with the infinite well conditions, ijj(x = 0) = ip(x = L) = 0. This model of an 
infinite well with a sinusoidal potential could simulate a quantum wire where 
the potential represents the atomic lattice within the wire. Thus, modelling 
the time evolution of a wave packet trapped within such a structure is a very 
simple "physical" test. 

5.1 Construction of the model 

In order to incorporate this potential into our previous Crank-Nicolson/finite 
element method we could simply calculate the average value of Eqn. (5.1) 
within each element and then use these for V e in the element equation given 
in Eqn. (4.22). In this way we obtain a constant FE approximation for the 
potential. However, a more general method would be to recalculate the po- 
tential element term in Eqn. (4.22) using the potential given in Eqn. (5.1). 
To obtain the best results we would need to use quadratic or higher order 
basis functions to take full advantage of this method of modelling the effects 
of the potential. However, for simplicity we will continue using linear ba- 
sis functions. Therefore, using this general method and linear elements, as 
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before, the potential element equation can be written as 



i r 

— -z / V(x)(j)ipdx 
* Jo 
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5> 
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~8h 



cos 



(1-0(1-0 (1-0(1 + 

(1 + 0(1-0 (i + oa + o 



where x ei is the position of node i of element e. Carrying out the 
and simplifying we obtain the potential element matrix 

iAL 8 



St 



8h k 3 (x ei 



x, 



Z2) 



C21 C22 



+1 
2 

(5.2) 
integral 

(5.3) 



where 

C n 

C\2 C21 
C22 



(2 — k 2 (x ei — x e2 ) 2 ) sin [kx ei ] — 2(sin [kx e2 ] + k(x ei — x e2 ) cos [kx ei ]) 
2 sin [kx €2 ] — 2 sin [fcx ei ] + k(x ei — x e2 )(cos [fcx ei ] + cos [kx e2 ]) 
(—2 + k 2 (x ei — x e2 ) 2 ) sin [fcx e2 ] + 2(sin [kx ei ] — k(x ei — x e2 ) cos [kx e2 ]). 

(5.4) 

also, if l e = x ei — x e2 we can write the full Schrodinger element equation 
(4.23) as 



"21" 


7 • f 6 ^ 


1 -1 " 




" C n 


Cl2 




1 2 


-1 1 


^21 


^22 


I 



^. (5.5) 

^ _6 e L W21 w ^ J ) 

In operator form this expression is given as 

A^ + i jB + cj ^ = 0, (5.6) 

with C representing the lattice-potential element matrix. Now, following the 
steps as in Section 4.1.2 we can write the Crank-Nicolson approximation as 

(A + i |B + c}) ij n+1 = (A - i |B + c}) (5.7) 
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Putting this in real form we have 

(A 7 + aB' + /?C) *p n+1 = (A' - aB f - (3C f ) (5.8) 

where A', B', a and (3 are as in Eqn. (4.28), but the potential term is 

-C n -C 12 

q, _ Cn C12 

—C21 —C22 

C21 C22 

and = iS%- 

Wave-packet with fc 7^ 

When we use a wave-packet with k = 2, 1000 spatial elements and dt = 0.1, 
combined with the potential parameters: A = 0.5 and fc = we have a 
packet that begins to move to the right with a small amount of "diffusion" to 
the left, Fig 5.2. As the packet interacts with many potential barriers in the 
lattice the probability conservation fluctuates and drops significantly after a 
very short time. Using a timestep size dt — 0.1 the probability, after 100 
timesteps, falls by 33% for 500 elements, 20% for 1000 elements and 12% for 
2000 elements, Fig. 5.3. From Fig. 5.4 it can be seen that as we increase the 
number of elements the damping effect diminishes, and for larger number of 
elements the timestep size has little effect. In Fig. 5.5 it can be seen that 
just reducing the timestep size has a worse effect, i.e. the errors due to a lack 
of resolving of the wave function at lower times are perpetuated along the 
time evolution, which results in highly unphysical results (probability > 1). 
Thus to obtain the best result we would need to use a time step of dt = 0.05 
(as reducing the timestep further would have little accuracy effect but would 
add greatly to the computation time) and use > 3000 spatial elements (more 



(5.9) 
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Figure 5.2: Wave packet with k = 2 spreads 
Im 2 are shown). 



over the lattice-potential (Re 2 + 
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Figure 5.3: Conservation of probability for various number of elements with 
dt — 0.1 (for packet with k = 2. Dashed line represent same results but for 
dt = 0.05 



elements would add to the resolving power for complicated wave functions 
at later times). 
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Figure 5.4: Probability at t = 100 against number of elements. Dashed line 
shows possible extrapolation, (for packet with k = 2) 
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Figure 5.5: Probability conservation for various timesteps, dt (for packet with 
k = 2). 
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Wave-packet with k = 

We also conducted the lattice potential simulation using a packet with k = 0, 
The simulation results, for 1000 elements and dt = 0.1, are shown in Fig. 5.6. 
It can be seen that the particle probability at the center of the well diffuses 
outwards from the central potential peak. However, as the packet starts with 
k = its energy is too low for it to tunnel through the lattice potential, so the 
probability of finding it near the well edges is very small even after t = 100 
timesteps. We can also note that the packet is more likely to be found within 
the central troughs - either side of the central peak of the lattice-potential, 
which is due to the packet being attracted to the states of lowest potential 
energy it is easily able to access. In this way we can confine the wave-packet 
for a long period within a certain location. In Fig. 5.7 we have a plot of 
probability conservation. This, again, shows that probability conservation 
falls slightly after t = 100 timesteps, which as before can be controlled by 
increasing the number of spatial elements. 
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Figure 5.6: Wave packet with k = is trapped within the central potential 
troughs (Re 2 + /m 2 are shown). 
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Figure 5.7: Conservation of area at each time step for packet with k G = in 
lattice-potential (for 1000 elements). 



Chapter 6 
Summary 



In Sec. 4.1 the Schrodinger equation for a wave packet in an infinite well was 
numerically solved by the use of linear finite elements for spatial variation 
and Crank-Nicolson for time evolution. Using a time step size dt = 0.5 we 
found that the probability was perfectly conserved even using as few as 50 
elements, Fig 4.3. However, it was seen that when a finite barrier is intro- 
duced within the well the conservation of probability is slightly "disturbed". 
When the wave packet interacts with the barrier, and reflected and trans- 
mitted waves are introduced, the probability conservation shows fluctuations, 
Figs. 4.9 and 4.10. These fluctuations are significantly reduced by increasing 
the number of elements, however reducing the timestep size below dt = 0.1 
has no effect. Thus we can conclude that these fluctuations are a result of 
the greater amount of detail in the wave function as it interacts with the bar- 
rier, which then requires more elements to resolve its shape. We also found 
that when a very high energy wave packet is used it can travel inside the 
well without noticing the barrier, hence no significant disturbances occur in 
its wave function, Fig. 4.13. The same is true for a very low energy wave 
function, i.e. it is trapped between the infinite well wall and the finite barrier 
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and so very little is transmitted, hence the low energy packet behaves as if 
it were trapped in a smaller infinite well. 

In Sec. 4.2.1 the Schrodinger equation, for a wave packet in an infinite 
well, was solved by applying linear continuous finite elements in space and 
time. In this method we obtain two non-unitary time stepping equations, 
Eqns. (4.34). One is explicit, as it is weighted with the two basis function 
components at time step t n , while the other is implicit, as it is weighted 
with the remaining basis function components at t n +\. As is expected the 
explicit case is unstable and it blows up, while the implicit case is heavily 
damped, Figs. 4. 14. However, we find that when we combine the explicit 
and implicit numerical equation we obtain the previous, Crank-Nicolson and 
finite element, time stepping equation. 

In Sec. 4.2.2 the Schrodinger equation is solved by applying the linear 
discontinuous space-time finite element method. Here we apply continuous 
finite elements in space and discontinuous finite elements in time. We again 
obtain a numerical equation that has a non-unitary form, Eqn. (4.44). We 
find that although the numerical solution obtained from this method suffers 
from damping it can be significantly controlled by decreasing the time step 
dt, and to a lesser extent by increasing the number of elements, Figs. 4.17 - 
4.19. From Fig. 4.20 it can be seen that by using 100 elements and a time 
step of dt = 0.01 we have almost zero damping, however this comes at the 
cost of greater computation time. 

In the case of introducing a finite potential barrier we find that the dis- 
continuous space-time method slightly reduces the size of fluctuations, which 
like the Crank-Nicolson method can be controlled by the number of elements 
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used, Fig. 4.21. The main issue with the discontinuous space-time method 
is still the damping. Although we can control this by reducing the timestep 
size, Fig. 4.22, this however significantly increases the computational cost. 
In Fig. 4.23 we have a direct comparison between the Crank-Nicolson and 
the discontinuous space-time methods for identical parameters: dt = 0.05 
and 250 elements. Even though the space-time method has less fluctuations 
it suffers from damping which causes a 1% drop in probability after 100 
timesteps. Another important factor is that the Crank-Nicolson method was 
five times faster in this test. 

As the Crank-Nicolson method came out on top we went on to investigate 
the limits of this method: i.e. the extent of the fluctuations which occur when 
the wave function interacts with a barrier/change of potential in a "physical" 
example. For this we modeled a particle in a quantum wire by using a wave 
packet in an infinite well containing a sinusoidal potential. From this we find 
that as the wave packet constantly interacts with the sinusoidal potential 
a loss of probably occurs i.e. a type of "damping" effect, Fig. 5.3. As for 
the simple barrier case, this is found to be an issue of resolving the wave 
function. After many timesteps the packet goes through a lot of transmissions 
and reflections, which in turn produces a complicated wave function. This 
complicated wave function then requires a greater number of elements to 
resolve the solution accurately. 

After the comparison of the novel space-time finite element method and 
the Crank-Nicolson method for the numerical solution of the Schrodinger 
equation we can conclude that in terms of accuracy, efficiency and conve- 
nience the Crank-Nicolson method is by far superior. 
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The logical extension to this work would be to go on and start to model 
a physical quantum device. However, this will need to be done in discrete 
steps, as a foundation of numerical tools has to be developed. The first step 
is to implement higher order spatial finite elements to model a wave function 
with greater accuracy, which could lead to a reduction of the fluctuations 
when the packet interacts with a change of potential. The next step would 
be to model more than one particle, and even introduce the interaction with 
electromagnetic fields. We will also need to implement different types of 
boundary conditions, for example periodic, open, absorbing, etc. Finally 
this work would need to be extended to two dimensions. Once these tools 
are developed we can apply them to model complicated quantum devices. 



Appendix A 

Quantum Physics BackGround 



A.l Wave-Particle Duality 

Einstein showed that the momentum of a photon is: 

P=~y (A.1) 

where p is the momentum, h is the Planck's constant, and A is the wavelength. 
This can be shown as follows. The energy of a photon is given as E — hv 
and its velocity as c = \v, where v is the photon frequency. Combining these 
we obtain: 

E = y (A.2) 

Now using Einstein's mass energy relation from the theory of relativity, E = 
mc 2 , we have: 

A = — , (A.3) 

mc 

where m is the relativistic mass of the photon 1 . Now noting that mass m 
multiplied by velocity c is momentum p, Eqn. (A.3) then becomes Eqn. (A.l). 
Therefore, this way it can be shown that electromagnetic wave packets of 



1 The rest mass of the photon is zero 
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evergy hv have a particle like behaviour. Since light can behave both as 
a wave, de Broglie reasoned in 1924 that matter also can exhibit this wave- 
particle duality. He further reasoned that matter would also obey Eqn. (A.l). 
In 1927, Davisson and Germer observed diffraction patterns by bombarding 
metals with electrons, confirming de Broglie's proposition. 

A. 2 Gaussian Wave Packets 

In terms of classical mechanics a particle's position and momentum can 
be measured exactly, however, in quantum mechanics this is not the case. 
The Heisenberg uncertainty principle implies that the highest precision with 
which the position and momentum of a particle can be measured is given by 
the minimum uncertainty relation: 

AxAP x = ^ (A.4) 

The wave function that satisfies this uncertainty relation is a Gaussian wave 
packet, 

ip(x) = (^2) 4 e lk " x e^ x - x ^' A(J \ (A.5) 

where x denotes the center of the wave packet, hk is the mean momentum 
of the packet, and a is the uncertainty in the position of the particle (Ax). 

In order to construct a numerical model of the wave packet we can sepa- 
rate the real and complex parts of the wave function: 

Mx) = I — ~ — 7z j 4 {cos(fcox) + i sm(k Q x)} e -(*-*o) 2 /4* 2 

This can now be written in the form of a vector equation, 

$(x) = 



Re[^(x)} 
Im[i/j(x)] 



cos(A;ox) 
sin(/c x) 



(A.6) 
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Figure A.l: Real and Imaginary parts of ^, with a = 2, k = 2, and x = 0. 

This vector notation of the complex wavefunction is the most convenient 
form for numerical computation. In Fig. A.l we have shown the real and 
imaginary plots of ij). 
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C++ Code Samples 



B.l Vector/Matrix Manipulation and the Vec- 
tor Library 

In this work we created a sub-program using the C++ vector class library 
in order to handle the large matrices and nodal vectors of the finite-element 
equations. One of the key advantages of vectors over arrays is their dynamic 
nature - once a vector is defined we can adjust its size at any time in or- 
der to handle differing amounts of data. Also writing a sub-program that 
constructs, handles and manipulates vectors would be a versatile tool which 
could be used for a variety of tasks. 

The vector sub-program was created as a header file (vector. h), and the 
first thing we had to do was implement the vector library and define a vector 
and matrix type: 

#include <vector> 

using namespace std; 

typedef vector<double> vectors; 

typedef vector<vectors> matrix; 

typedef vectors :: iterator vecit; 

Once the vector and matrix types were defined we created a vector (vec) and 
Matrix (mat) classes, which incorporated member function (or methods) to 
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take care of some important vector and matrix tasks: 

class vec:public vectors{ //Vector Class 
public : 

vecQO //General vector constructor 

vec(int i){ 

for(int k=0;k<i ;k++){ //Constructor for vector with all zero components 

push_back(0) ; 

} 

} 

void print (); //Function to output vector 

void pop_front(); //Removes the first element of a vector (pop_back() exists in standard library) 

void pop_ends(); //Removes the first and last component (required for boundary conditions) 

}; 

void vec : :print () { 

int i=size () ; 

cout<< M \n\n M ; 

for(int k=0;k<i;k++){ 

cout«*(begin()+k)<<" "; 

} 

cout<< M \n\n M ; 
} 

void vec : :pop_f ront () { 
int i=size()-l; 
for(int k=0;k<i;k++){ 
*(begin()+k)=*(begin()+k+l) ; 
} 

pop_back() ; 
} 

void vec: :pop_ends(){ 
pop_f ront () ; 
pop_back() ; 
} 

class mat: public matrix{ //Matrix Class 
public : 

matQO //General matrix constructor 

mat (int i, int j){ //Constructor for matrix with all zero components 
vec row; 

for(int k=0;k<j ;k++){ 
row.push_back(0) ; 
} 

for(int l=0;Ki;l++){ 

push_back(row) ; 

} 

} 

void print (); //Function to output vector 
unsigned int rows(); //Outputs number of rows in matrix 
unsigned int cols(); //Outputs number of columns in matrix 
friend void pop_front(); 

void pop_edges(); //Removes the outer components of matrix (required for boundary conditions) 
}; 

unsigned int mat:: rows (){ 

return size () ; 

} 

unsigned int mat::cols(){ 
return (* (begin() ) ) . size () ; 
} 

void mat :: print (){ 
int i=rows () ; 
int j=cols () ; 
vec row( j ) ; 
cout<< M \n\n M ; 
for(int k=0;k<i;k++){ 

for(int l=0;Kj ; 1++) cout<< (* (begin() +k) ) [1]« M "; 

cout<< M \n\n M ; 

} 

} 

void mat : :pop_edges(){ 
int i=size()-l; 
for(int k=0;k<i;k++){ 
*(begin()+k)=*(begin()+k+l) ; 
} 

pop_back() ; 
pop_back() ; 
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for(int j=0; j<i-l; j++){ 
for (int k=0;k<i;k++){ 

*((*(begin()+j)) .begin()+k)=* ( (* (begin()+j ) ) .begin()+k+l) ; 
} 

(*(begin()+j)) .pop_back() ; 
(*(begin()+j)) .pop_back() ; 
} 
} 



In the vector. h library we also included some important functions and over- 
loaded operators: 

//Dot product of two vectors x and y of equal size and return a double value 
double operator* (vec &x, vec &y){ 
if (x. size () ! =y . size () ){ 

cout<<x. size () «"\n" ; 
cout<<y . size () «"\n" ; 

cout<< M \n\n Vectors of unequal size M <<x . size () <<" and "<<y . size () «"\n" ; 
return 0; 

} 

else{ 

int max=x. size ()-l; 
double dot=0; 
for(int i=0;i<=max;i++){ 
dot+=x[i]*y[i] ; 

} 

return dot ; 

} 

} 

//Matrix multiplication of two matrices a and b of sized ixj and jxk and return a matrix of size ixk 
mat operator* (mat &a, mat &b){ 
if (a.colsO !=b.rows()){ 

cout<< M \n\n M << M Matrix sizes not compatible: "«a.rows()«" x "<< a.colsO 

«" and "«b.rows()<<" x M «b . cols () «"\n\n" ; 
mat x ( 1 , 1 ) ; 
return x; 

} 

else{ 

int maxi=a.rows () ; 
int maxk=b . cols () ; 
int maxj=a. cols () ; 
vec row; 
mat c; 
c . clear () ; 
double element; 
for(int i=0;i<maxi;i++){ 
row. clear () ; 

for (int k=0;k<maxk;k++){ 
element=0; 

for(int j=0; j<maxj ; j++)element+=(a[i] [j]*b[j] [k] ) ; 
row .push_back(element) ; 

} 

c .push_back(row) ; 
} 

return c; 
} 

} 

//Vector addition of two vectors x and y of size i 
vec operator+ (vec &x, vec &y){ 
if (x. size () ! =y . size () ){ 

cout<< M \n\n M << M Vectors sizes not compatible: "<<x.size() 

<<" and "<<y.size(); 
vec z(l) ; 
return z; 

} 

else{ 

vec z; 

int maxi=x . size () ; 
z . clear () ; 

for(int i=0;i<maxi;i++){ 

z . push_back (x [i] +y [i] ) ; 
} 

return z; 
} 

} 
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//Matrix addition of two matrices a and b of size ixj 
mat operator+ (mat &a, mat &b){ 

if (a. cols () ! =b. cols ()&&a. rows () !=b.rows()){ 

cout<< M \n\n M << M Matrix sizes not compatible: "«a.rows()«" x "<< a.colsQ 

«" and M «b.rows()<<" x M «b.cols(); 
mat c ( 1 , 1 ) ; 
return c; 

} 

else{ 

mat c; 
vec row; 

int maxi=a.rows () ; 
int maxj=a. cols () ; 
c . clear () ; 

for(int i=0;i<maxi;i++){ 
row. clear () ; 

for(int j=0; j<maxj ; j++){ 

row .push_back(a[i] [j] +b [i] [j] ) ; 

} 

c.push_back(row) ; 

} 

return c; 

} 

} 

//Vector subtraction of two vectors x and y of size i 
vec operator- (vec &x, vec &y){ 
if (x. size () ! =y . size () ){ 

cout<< M \n\n M << M Vectors sizes not compatible: "<<x.size() 

<<" and "<<y.size(); 
vec z(l) ; 
return z; 

} 

else{ 

vec z; 

int maxi=x . size () ; 
z . clear () ; 

for(int i=0;i<maxi;i++){ 

z . push_back (x [i] -y [i] ) ; 
} 

return z; 
} 

} 

//Matrix subtraction of two matrices a and b of size ixj 
mat operator- (mat Sea, mat &b){ 

if (a. cols () ! =b . cols ()&&a. rows () ! =b .rows () ) { 

cout<< M \n\n M << M Matrix sizes not compatible: "«a.rows()«" x "<< a.colsQ 

«" and M «b.rows()<<" x M «b.cols(); 
mat c ( 1 , 1 ) ; 
return c; 

} 

else{ 

mat c; 
vec row; 

int maxi=a.rows () ; 
int maxj=a. cols () ; 
c . clear () ; 

for(int i=0;i<maxi;i++){ 
row. clear () ; 

for(int j=0; j<maxj ; j++){ 

row .push_back(a[i] [j] -b [i] [j] ) ; 

} 

c.push_back(row) ; 

} 

return c; 

} 

} 

//Matrix transpose 

mat transpose (mat a){ 

int i=a.rows () ; 
int j=a. cols () ; 
int k, 1; 
mat t ( j , i) ; 
for(k=0;k<i;k++){ 

for(l=0;Kj ;l++){ 

t[l] [k]=a[k] [1]; 

} 

} 

return t ; 
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//Multiplication of a scalar with a matrix 
mat operator* (double &a, mat &b){ 
mat c; 



vec row; 

int maxi=b .rows () ; 
int maxj=b . cols () ; 
c . clear () ; 

for(int i=0;i<maxi;i++){ 
row. clear () ; 

for(int j=0; j<maxj ; j++){ 

row .push_back(a*b [i] [j] ) ; 

} 

c.push_back(row) ; 

} 

return c; 



//Matrix and vector multiplication 

vec operator* (mat &a,vec &v){ 
int max=v . size () ; 
vec x(max) ; 



for (int i=0; i<max; i++){ 

for(int j=0; j<max; j++){ 

x[i]+=a[i] [j]*v[j] ; 

} 

} 

return x; 



B.2 LU Decomposition and Forward and Back- 
ward Substitution 

In this work the global matrix equations, 

A.X = (L-U)-X = B, 

are solved for the nodal values, X, by using simple LU decomposition as 
described in [15]. The biggest limitation with this simple method is that the 
entire A matrix is used even though it is sparse diagonal and thus contains 
a large number of zeros. A decision was made to use this method as it was 
easy to program, and the main aim of the work was to compare the results 
rather than actually build a fast program. A library file, matrixop.h, was 
created to hold the decomposition and substitution routines so they could be 
easily installed into new simulations. The LU decomposition routine takes a 
matrix a[i][j] and replaces it by the LU decomposition of itself: 

//LU Decomposition 

void lu(mat &a, vec feindx, double &d){ 
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const double tiny=l . Oe-20 ; 

int i,imax, j ,k; 

double big,dum, sum, temp; 

int n=a. rows () ; 

vec vv(n) ; 

d=1.0; 

for(i=0;i<n;i++){ 
big=0.0; 

for(j=0; j<n; j++){ 

if ((temp=fabs(a[i] [j] ) ) >big)big=temp ; 
} 

if(big==0.0) cout<< M Singular matrix"; 

vv[i]=1.0/big; 

} 

for(j=0; j<n; j++){ 
for(i=0;i<j ;i++){ 
sum=a[i] [j] ; 

for(k=0;k<i;k++) sum -= a[i] [k] *a[k] [j] ; 

a[i] [j]=sum; 

} 

big=0.0; 

for(i=j ;i<n;i++){ 
sum=a[i] [j] ; 

for(k=0;k<j ;k++) sum -= a[i] [k] *a[k] [j] ; 
a[i] [j]=sum; 

if ( (dum=vv [i] *f abs (sum) ) >=big) { 

big=dum; 

imax=i ; 

} 

} 

if (j !=imax){ 
for(k=0;k<n;k++){ 
dum=a[imax] [k] ; 
a[imax] [k]=a[j][k]; 
a[j] [k]=dum; 
} 

d=-d; 

vv [imax] =vv [ j ] ; 
} 

indx [j]=imax; 

if (a[j] [j]==0.0) a[j] [j]=tiny; 
if (j!= n-l){ 
dum=1.0/(a[j][j]); 

for(i=j+l;i<n;i++) a[i] [j] *= dum; 

} 

} 

} 



The routine for forward and backward substitution solves the equation A • 
X = B. Here a[i][j] is input as the LU decomposed version of the matrix A, 
and b[j] is input as the right hand vector B. Then on output the results for 
X are returned in B: 



//LU Substitution 

void luback(mat &a, vec &indx,vec &b){ 
int i, ii=0,ip,j; 
double sum; 
int n=a. rows () ; 

for(i=0;i<n;i++){ 
ip=indx[i] ; 
sum=b [ip] ; 
b[ip]=b[i] ; 
if (ii!=0){ 

for(j=ii-l; j<i; j++) sum -= a [i] [j] *b [j] ; 
} 

else if(sum != 0.0){ 

ii = i + 1; 

} 

b [i] =sum; 
} 
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for(i=n-l;i>=0;i--){ 
sum=b [i] ; 

for(j=i+l; j<n; j++) sum -= a[i] [j] *b[j] ; 

b [i] =sum/a[i] [i] ; 

} 

} 



Bibliography 



[1] WIAS, Quantum Mechanical and Macroscopic Models for Optoelec- 
tronic Devices, 

http : //www . wias-berlin . de/pro j ect-areas/micro-el/f zt/FZT86-D4 . html 

[2] R. L. Liboff, Introductory Quantum Mechanics, Addison- Wesley Pub- 
lishing Company (1993). 

[3] Wikipedia, Quantum Mechanics, 

http : //en. wikipedia . org/ w/ index . php?t itle=Quantum_mechanics 
&oldid=52542410 

[4] R. Chen, Z. Xu and L. Sun, Finite- Difference Scheme to Solve 
Schrddinger Equations, Phys. Rev. E 47, 3799-3802 (1993). 

[5] Wikipedia, Finite Element Method, 

http : //en . wikipedia . org/wiki/Finite_element_method 

[6] O. C. Zienkiewicz and Y. K. Cheung, The Finite Element Method in 
Engineering Science, McGraw-Hill London (1967). 

[7] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Ele- 
ment Methods, Springer- Verlag (2002). 

99 



BIBLIOGRAPHY 



100 



[8] W. B. Bickford, Finite Element Method, Irwin Publishing (1994). 

[9] C. F. Gerald and P. O. Wheatley, Applied Numerical Analysis, Pearson 
Addison Wesley (2003). 

[10] W. Dettmer and D. Peric, Analysis of time integration algorithms for 
the finite element solutions of incompressible Navier-Stokes equations, 
Comp. Meth. App. Mech. 192 1177-1226 (2003). 

[11] N. Watanabe and M. Tsukada, Finite element approach for simulating 
quantum electron dynamics in a magnetic field, Jour. Phys Soc. Jap. 69 
No.9 2962 (2000). 

[12] L. Ramda Ram-Mohan, Finite Element and Boundary Element Applica- 
tions in Quantum Mechanics, Oxford Texts in Applied and Engineering 
Mathematics (2002). 

[13] L. Solymar and D. Walsh, Electrical Properties of Materials, Oxford 
University Press (2004). 

[14] IBM, Scanning Tunnelling Microscopy, 

http : / / www . almaden . ibm . com/ vis/ stm/ stm . html 

[15] William H. Press, Numerical Recipes in C++: The Art of Scientific 
Computing , Cambridge University Press (2002). 



List of Figures 



2.1 Infinite Square well 12 

2.2 Particle in a box wavefunctions 13 

2.3 Hermite polynomials and the corresponding energy eigenstates. 14 

2.4 Quantum harmonic oscillator eigenstates 14 

3.1 Node numbering and coordinates of rectangular element. ... 32 

3.2 Linear discontinuous finite elements in space-time 37 

4.1 Wave packet trapped in a infinite well (cont. on next page). . 46 

4.2 (Cont. from previous page) Wave packet trapped in an infinite 

well 47 

4.3 Conservation of area at each time step for varying number of 
elements 48 

4.4 Wave packet with k = spreads over the infinite well (Re 2 + 

Im 2 are shown) 49 

4.5 Conservation of area at each time step for packet with k = 0. 50 

4.6 Infinite potential well with a barrier 50 

4.7 Wave packet trapped in a infinite well with barrier at elements 

120 < x e < 130 (cont. on next page) 54 

4.8 (Cont. from previous page) Wave packet trapped in an infinite 

well with barrier at elements 120 < x e < 130 55 



101 



LIST OF FIGURES 



102 



4.9 Conservation of area at each time step for varying number of 
elements. (Also on the same plot is area conservation for 500 
elements and time step dt=0.1) 56 

4.10 Conservation of area at each time step for varying size of dt. 
(Also on the same plot is area conservation for 500 elements 

and time step dt=0.1) 57 

4.11 Wave packet with initial wave vector k = 3 trapped in a 
infinite well with barrier at elements 120 < x e < 130 59 

4.12 Wave packet with initial wave vector k = 1 trapped in a 
infinite well with barrier at elements 120 < x e < 130 60 

4.13 Conservation of area at each time step for varying initial wave 
vector k 61 

4.14 Area after each time-step for explicit scheme (left) and im- 
plicit scheme (right). Area explodes for the explicit case, and 
remains finite but not constant for the implicit case 63 

4.15 Conservation of probability at each time step for 100 elements 

and varying size of timestep size: dt 67 

4.16 Conservation of probability at each time step for dt = 0.05 

and varying number of elements 68 

4.17 Conservation of probability at t = 100 against dt for 100 ele- 
ments 69 

4.18 Close-up of conservation of probability at t = 100 against dt 

for 100 elements 69 

4.19 Conservation of probability at t = 100 against number of ele- 
ments for dt = 0.05 70 

4.20 Conservation of probability at each time step for 50 and 100 
elements and timesteps dt = 0.1, dt = 0.05, and dt = 0.01. . . 70 



LIST OF FIGURES 103 

4.21 Conservation of probability at each timestep, for various num- 
ber of elements (timestep size dt = 0.01) 72 

4.22 Conservation of probability at each timestep, for various timestep 
sizes dt (100 elements) 73 

4.23 Conservation of probability at each timestep, for discontinuous 
space-time (DST) and Crank-Nicolson (CN) methods for 250 
elements and timestep dt = 0.05 74 

5.1 Model of a periodic lattice potential within a quantum wire. . 76 

5.2 Wave packet with k = 2 spreads over the lattice-potential 
(Re 2 + Im 2 are shown) 79 

5.3 Conservation of probability for various number of elements 
with dt = 0.1 (for packet with k = 2. Dashed line represent 
same results but for dt = 0.05 80 

5.4 Probability at t = 100 against number of elements. Dashed 

line shows possible extrapolation, (for packet with k = 2) . . 81 

5.5 Probability conservation for various timesteps, dt (for packet 
with k = 2) 81 

5.6 Wave packet with k = is trapped within the central poten- 
tial troughs (Re 2 + Im 2 are shown) 83 

5.7 Conservation of area at each time step for packet with k = 

in lattice-potential (for 1000 elements) 84 

A.l Real and Imaginary parts of with a = 2, k = 2, and x = 0. 91 



